function Nm=newton_bb(dofs,mesh,f)
% f(x,u) is nonlinear function of newton item
V=mesh.V;
T=mesh.T;
d=mesh.d;
posT=size(T,1);

pos=1;

for k = 1:posT
    V1=V(T(k,1),:);V2=V(T(k,2),:);V3=V(T(k,3),:);
    x13 = V1(1)-V3(1);
    y13 = V1(2)-V3(2);
    x23 = V2(1)-V3(1);
    y23 = V2(2)-V3(2);
    J = x13*y23 - x23*y13; 
    [qp,qw,b1,b2,b3]=quadrule;
     qp_glob = zeros(size(qp));
%       qp_glob(:,1) = V3(1) + x13*qp(:,1) + x23*qp(:,2);
%       qp_glob(:,2) = V3(2) + y13*qp(:,1) + y23*qp(:,2);  
      % evaluate at integrate points, 
        qp_glob(:,1) = b1*V1(1) + b2*V2(1) + b3*V3(1);
      qp_glob(:,2) =  b1*V1(2) + b2*V2(2) + b3*V3(2); 
         loc_dof = get_tri_dof(dofs,k);
         c_loc=mesh.c(loc_dof);
   u_loc=vdm23(d,b1,b2,b3)*c_loc;
   um=feval(f,u_loc);
   um=diag(um.*qw);
   nm=vdm23(d,b1,b2,b3)'*um*vdm23(d,b1,b2,b3);
   [i,j,s]=find(nm);
   L = length(i);
    Indx1(pos:(pos + L-1)) = dofs(k,i);
    Indx2(pos:(pos + L-1)) = dofs(k,j);
    S(pos:(pos + L-1)) = J*s/2;
    pos = pos + L;
   end
dim = max(max(dofs));
Nm = sparse(Indx1(1:(pos-1)),Indx2(1:(pos-1)),S(1:(pos-1)),dim,dim);